#########
#### Gay rights analysis for Warshaw CUP chapter
#### Note: The first third of code builds MRP estimates using LMER in R. The point estimates 
#### from this model are accurate, but the uncertainty is too small. So in the next third of 
#### the code I run a "production" version MRP model using Stan.  I also check that both the 
#### LMER and Stan-based point estimates are identical.  Finally, the last third of the 
#### code runs generates analyses for the CUP Data Science chapter text.
#########
library(car)
library(plyr); library(dplyr)
library(foreign)
library(ggplot2)
library(lme4)
library(mapproj)
library(maps)
library(parallel)
library(pbapply)
library(pscl)
library(reshape2)
library(rms)
library(arm)
library(rstan)
library(stringr)


#########
#### Build graphs (Figures 1 and 2) and analyses for text
#########

cities<-read.csv(file="tw_apsr_cities_finaldata.csv",colClasses = "character")
gayrights<-read.csv(file="cities_gayrights_policies.csv",colClasses = "character")
gayrights_po<-read.csv(file="cities_posteriors_gaymarriage.csv",colClasses = "character")
colnames(gayrights_po)[which(colnames(gayrights_po)=="mrp_stan")]<-"gayrights.stan"
colnames(gayrights_po)[which(colnames(gayrights_po)=="mrp_lmer")]<-"gayrights.lmer"
colnames(gayrights_po)[which(colnames(gayrights_po)=="raw")]<-"gayrights.raw"
states<-read.csv(file="states_icpsr.csv",colClasses = "character")
states<-states[,c(1,2)]
states$NAME <- gsub(" ","",states$NAME)
states$NAME[states$NAME=="DistrictofColumbia"] <- "District of Columbia"
states$NAME[states$NAME=="NewHampshire"] <- "New Hampshire"
states$NAME[states$NAME=="NewJersey"] <- "New Jersey"
states$NAME[states$NAME=="NewMexico"] <- "New Mexico"
states$NAME[states$NAME=="NewYork"] <- "New York"
states$NAME[states$NAME=="NorthCarolina"] <- "North Carolina"
states$NAME[states$NAME=="NorthDakota"] <- "North Dakota"
states$NAME[states$NAME=="RhodeIsland"] <- "Rhode Island"
states$NAME[states$NAME=="SouthCarolina"] <- "South Carolina"
states$NAME[states$NAME=="SouthDakota"] <- "South Dakota"
states$NAME[states$NAME=="WestVirginia"] <- "West Virginia"


cities<-merge(cities,states,by.x="abb",by.y="STATEAB")
cities2<-merge(cities,gayrights,by.x=c("NAME", "city"),by.y=c("State","City"), all.x=T)
cities2<-merge(cities2,gayrights_po,by.x=c("city_id"),by.y=c("city_id"), all.x=T)
cities2$Gay.Employer.Score<-as.numeric(as.vector(cities2$Gay.Employer.Score))
cities2$policy<-as.numeric(as.vector(cities2$policy))
cities2$gayrights.stan<-as.numeric(as.vector(cities2$gayrights.stan))
cities2$gayrights.lmer<-as.numeric(as.vector(cities2$gayrights.lmer))
cities2$gayrights.raw<-as.numeric(as.vector(cities2$gayrights.raw))
cities2$pres_2008<-as.numeric(as.vector(cities2$pres_2008))
cities2$stan.lower<-as.numeric(as.vector(cities2$stan.lower))
cities2$stan.upper<-as.numeric(as.vector(cities2$stan.upper))
cities2$mrp_ideology<-as.numeric(as.vector(cities2$mrp_ideology))

cities2<-subset(cities2, city_pop>.5)

###########
## Figure 2: Association between public opinion and gay rights policies
#############


cities2<-subset(cities2, !is.na(Gay.Employer.Score))
cities2<-cities2[order(cities2$gayrights.stan),]
cities2$hj<-0
cities2$hj[cities2$city_pop>5] <- rep(c(-1.65,1.65),
                                      length.out=dim(cities2[cities2$city_pop>5,])[1])

cities2$city2<-cities2$city
cities2$city_pop<-as.numeric(as.vector(cities2$city_pop))
cities2$city2[cities2$city_pop<6]<-""


library(ggplot2)
c <- ggplot(cities2, aes( gayrights.stan,Gay.Employer.Score, label=city2))
c <- c + geom_point(size=.5)
c <- c + geom_text(size=3,aes(vjust=hj))
c <- c + stat_smooth(method = "loess", formula = y ~ x, size = 1)
c <- c + scale_y_continuous(name="Index of Gay Rights for Municipal Employees")
c <- c + scale_x_continuous(name="Support for Same-Sex Marriage",breaks=c(.4,.5,.6,.7,.8,.9), 
labels=c("40%", "50%", "60%", "70%", "80%", "90%"),lim=c(.35,.85))
c <- c + theme(axis.title.x = element_text(face="bold",  size=14))
c <- c + theme(axis.title.y = element_text(face="bold",  size=14))
c <- c + theme_bw()
pdf(file="figures/figure3_cities_gayrights_responsiveness.pdf", width=8, height=10)
c
dev.off()

